######################################################################
# Time-stamp: <05/09/09 11:06:21 ostolop>
######################################################################
package EP::EPC::Core::GO_Annotation;

use strict;

use EP::Config;
use EP::Common;
use EP::Common::DB;
use EP::Common::Error qw ( :try );
use EP::Common::General;

use EP::EPC::Core::Component;

use CGI qw/:standard/;

use GO::TermFinder;
use GO::AnnotationProvider::AnnotationParser;
use GO::OntologyProvider::OntologyParser;
use GO::View;

use GO::Utils::General qw (CategorizeGenes);

use File::Temp;
use File::Basename;

our @ISA = qw ( EP::EPC::Core::Component );

$GO::TermFinder::WARNINGS = 0;

# sub: process_component
#   Processes the GO_Annotation component
#   It expects the query parameters gene_ids (a list of gene names/ids), multtest_correction and cutoff to be set.
#
# returns: 
#   OK if everything went fine
sub process_component {
  my $self = shift;

  my @gene_ids_to_map;

  my $gene_ids = $self->{query}->param('gene_ids');
  if ( $gene_ids ) {
    $gene_ids =~ s/^\s+|\s+$//g;	# get rid of leading/trailing spaces
    @gene_ids_to_map = split /\n|\r\n|\r/, $gene_ids;	# split on newlines
  }

  my $selected_cluster = $self->{query}->param ( "current_cluster_selection" );
  if ( $selected_cluster ) {
    my $gene_ids_from_cluster = $self->get_gene_ids_from_cluster ( $selected_cluster );
    push @gene_ids_to_map, @{$self->{src_dataset}->data->map_gene_ids_for_GO ( $gene_ids_from_cluster ) };
  } else {
    # if no cluster selected, use the whole selected dataset
    push @gene_ids_to_map, @{$self->{src_dataset}->data->map_gene_ids_for_GO ( ) };
  }

  my $result = $self->map_genes_GO_terms ( \@gene_ids_to_map,
					   $self->{query}->param('multtest_correction'),
					   $self->{query}->param('cutoff') );

  my $dataset_folder_path = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $self->{src_dataset}->folder->location;
  my $output_www_base = "$EP::Config::h->{EPNG_USRDATA_URL}/" . $self->{src_dataset}->folder->location . "/";

  my $htmlfile = new File::Temp ( DIR => $dataset_folder_path, SUFFIX => ".html", UNLINK => 0 );
  my $htmlfile_url = $output_www_base . basename ( $htmlfile );

  print $htmlfile $result;
  close $htmlfile;

  my $process_result = { html => $htmlfile_url,
			 src_dataset_id => $self->{src_dataset}->id,
			 dst_dataset_id => $self->{src_dataset}->id }; # HACK: produce HTML output and place it all directly in <html></html> as a result!

  $self->{CD}->{process_result} = $process_result;

  return $self->record_analysis_history();
}

sub read_tree_data_from_file {
  my ($file) = shift;

  my @DATA;

  open ( CLUSTER, $file ) or die "Couldn't open the hierarchical clustering file $file: $!\n";
  while (<CLUSTER>) {
    next if /^\s*#/;

    # From the original clustering output
    my @data_row = split /\s+/;
    push @DATA, \@data_row;
  }
  close CLUSTER;

  return \@DATA;
}

# sub: initialize
#   Initialize the GO_Annotation object
#   The initialize routine will read in the ontology file belonging to the GO category and
#   will then read the annotation file that maps the genes of a certain organism to the ontology
#   tree.
#   The initialize routine expects the query paramters go_category and dataset_species to be set.
#
# returns:
#   Object of type GO_Annotation
#
# see also:
#   <EP::EPC::Core::Dispatcher::initialize>
sub initialize {
  my $self = shift;
  $self->SUPER::initialize(@_);

  throw EP::Common::Error ( -value => AFFY_TRANSFORM_RAW_DATA ) if $self->{src_dataset}->type eq "affy";

  select ( STDERR );

  $self->{category} = $self->{query}->param('go_category');
  $self->{species} = $self->{src_dataset}->species;
  $self->{aspect} = uc(substr($self->{category}, 0, 1));

  if ( !exists( $EP::Config::h->{SPECIES_METADATA}{$self->{species}}{annotationFile} ) ) {
      # No annotation file present for current dataset species
      throw EP::Common::Error ( -text => "No Gene Ontologies present for selected dataset, species " . $self->{species} );
  }

  # obtain the general ontology structure file
  $self->{session}->status ( "Reading Ontology Tree ..." );
  my $ontology_file = "$EP::Config::h->{EPNG_GODATA_PATH}/$self->{category}.ontology";
  if ( -f $ontology_file ) {
    if ( -f "$ontology_file.obj" &&
	 -M "$ontology_file.obj" < -M $ontology_file ) {
      # serialized objects should be there and newer than ontology file in case of updates
      warn "Reading from stored file" if DEBUG;
      $self->{ontologyProvider} = GO::OntologyProvider::OntologyParser->new ( objectFile => "$ontology_file.obj" );
    } else {
      warn "Serializing ontology to $ontology_file.obj..." if DEBUG;
      $self->{ontologyProvider} = GO::OntologyProvider::OntologyParser->new ( ontologyFile => $ontology_file );
      $self->{ontologyProvider}->serializeToDisk ( filename => "$ontology_file.obj" );
    }
  } else {
    die "Could not find ontology file for " . $self->{category};
  }

  # obtain the mapping file for this species
  $self->{session}->status ( "Reading Annotation Mappings ..." );
  my $annotation_file = $EP::Config::h->{EPNG_GODATA_PATH} . "/" . $EP::Config::h->{SPECIES_METADATA}->{$self->{species}}->{annotationFile};
  if ( -f $annotation_file ) {
      if ( -f "$annotation_file.obj" &&
	 -M "$annotation_file.obj" < -M $annotation_file ) {
	# serialized objects should be there and newer than annotation file in case of updates
      warn "Reading from stored file" if DEBUG;
      $self->{annotationProvider} = GO::AnnotationProvider::AnnotationParser->new ( objectFile => "$annotation_file.obj" );
    } else {
      warn "Serializing annotation to $annotation_file.obj..." if DEBUG;
      $self->{annotationProvider} = GO::AnnotationProvider::AnnotationParser->new ( annotationFile => $annotation_file );
      $self->{annotationProvider}->serializeToDisk ( filename => "$annotation_file.obj" );
    }
  } else {
    die "Could not find the annotation file for species: " . $self->{species};
  }

  # Create the background distribution set, the currently selected dataset
  # GO TermFinder wants completely unambiguous gene names for the population input
  my (@ambiguous, @population, @notFound);
  my @genes = @{ $self->{session}->current_root_dataset()->data()->map_gene_ids_for_GO };

  CategorizeGenes ( annotation  => $self->{annotationProvider},
		    genes       => \@genes,
		    ambiguous   => \@ambiguous,
		    unambiguous => \@population,
		    notFound    => \@notFound);

  throw EP::Common::Error ( -text => "Couldn't form a background set of genes for GO annotation: check that you have appropriate identifiers available." ) if scalar @population == 0;

  warn "Background initialization: " . scalar @genes . " genes total, " . scalar @ambiguous . " ambiguous, " . scalar @notFound . " not found.";

  $self->{termFinder} = GO::TermFinder->new( annotationProvider => $self->{annotationProvider},
					     ontologyProvider => $self->{ontologyProvider},
					     population => \@population,
					     aspect => $self->{aspect} );
}

# sub: map_genes_GO_terms
#   Maps a list of genes to the GO annotation if possible. If genes can not be mapped (due to incorrect/dubious gene names), they
#   will be put into seperate arrays and not used for finding significant GO terms. Those genes that can be mapped are then searched
#   for a significant overrepresentation of GO categories using a hypergeometric testing procedures. Multiple hypothesis correction can be
#   performed if neccessary.
#
# parameters:
#   genes - List of gene names/ids
#   multtest_correction - Type of multiple hypothesis correction, can be any of: bonferroni, simulate, none
#   cutoff - P-value significance cutoff
#
# returns:
#   Html formatted results. 
#   Both a graphical representation of the GO tree and a table with the significant GO categories above the cutoff are shown.
sub map_genes_GO_terms {
  my $self = shift;
  my ($genes, $multtest_correction, $cutoff) = @_;
  my (@list, @notFound, @ambiguous);

  select (STDERR);

  $self->{session}->status ( "Categorizing genes ..." );

  CategorizeGenes ( annotation  => $self->{annotationProvider},
		    genes       => $genes,
		    ambiguous   => \@ambiguous,
		    unambiguous => \@list,
		    notFound    => \@notFound);

  # no matching gene names found
  if (!@list) {
    return "<body><p><i>None of the gene names were recognized</i></p></body>";
  }

  $self->{session}->status ( "Finding significant GO terms ..." );
  my @pvalues = $self->{termFinder}->findTerms ( genes        => \@list,
						 correction   => $multtest_correction,
						 calculateFDR => 1);	

  # get the number of genes above the cutoff
  my $categories_above_cutoff = 0;
  foreach my $pvalue (@pvalues) {
      if ( $multtest_correction eq 'none' ) {
	  # if no multiple testing correction is performed, make sure
	  # we set the corrected pvalue to the original, otherwise GO::View and others go wrong
	  $pvalue->{CORRECTED_PVALUE} = $pvalue->{PVALUE};
      }
      $categories_above_cutoff++ if $pvalue->{CORRECTED_PVALUE} <= $cutoff;
  }

  my $HTML_OUTPUT = start_html( -title => "GO",
				-target => "_blank",
				-bgcolor => "white"
			      );

  $HTML_OUTPUT .= h2( {-align => "center"}, 
		      "GO TermFinder results");


  if ( $categories_above_cutoff == 0 ) {
    $HTML_OUTPUT .= h3( { -align => "center" }, "No significant GO categories were found above $cutoff cutoff for the " . scalar @list . " recognized genes." );

    # unused gene names
    if (@ambiguous) {
      $HTML_OUTPUT .= "<p>The following gene(s) are ambiguously named, and so were not used:<br>\n";
      $HTML_OUTPUT .= join(" ", @ambiguous) . "</p>\n";
    }
    if (@notFound) {
      $HTML_OUTPUT .= "<p>The following " . scalar @notFound . " gene(s) were not recognized, and so were not considered:<br>\n";
      $HTML_OUTPUT .= "<pre>" . join("<br>", @notFound) . "</p></pre>\n";
    }

    $HTML_OUTPUT .= end_html();

    return $HTML_OUTPUT;
  }

  # Some general output settings
  my $username = $self->{current_user}->username;
  my $resultDir = $EP::Config::h->{EPNG_USRDATA_PATH} . "/" . $username;
  my $resultWWW = $EP::Config::h->{EPNG_USRDATA_URL} . "/" . $username;
  my $rand = "$$" . rand(100);
  my $imageFile = $self->{src_dataset}->filename . ".vis." . $rand;

  $self->{session}->status ( "Creating output ..." );

  # Generate the GO tree
  my $params = { -ontologyProvider   => $self->{ontologyProvider},
		 -annotationProvider => $self->{annotationProvider},
		 -termFinder         => \@pvalues,
		 -aspect             => $self->{aspect},
		 -configFile         => "$EP::Config::h->{CONF}/GO_TermFinder.conf",
		 -imageUrlRoot       => $resultWWW,
		 -imageDir           => $resultDir,
		 -imageName          => $imageFile,
		 -imageLabel         => "GO tree view Expression Profiler",
		 -nodeUrl            => $EP::Config::h->{SPECIES_METADATA}->{$self->{species}}->{nodeUrl},
		 -geneUrl            => $EP::Config::h->{SPECIES_METADATA}->{$self->{species}}->{geneUrl},
		 -pvalueCutOff       => $cutoff };

  my $goView = GO::View->new( %$params );

  my $vis_format = 1; # $self->{current_user}->visualization_format;


  # 0 - SVG, 1 - PNG
  if ($vis_format == 0) {
      $HTML_OUTPUT .= table ( {-border => 0, -align=>"center", -width=>"100%"},
			      caption(''),
			      Tr ( { -align => "CENTER", -valign => "TOP" },
				   [
				    td (
					{ -style => "border: thin solid black; padding: 2px" },
					div (
					     { -style => "overflow:auto; height:500px; width:100%" },
					 svg_html ( "$resultWWW/$imageFile.svg", "100%", "100%" )
					     )
					)
				    ]
				   )
			      );
  }

      if ($goView->graph) {
	  if ($vis_format == 0) {
	      $goView->graph->as_svg("$resultDir/$imageFile.svg");
	  } else {
#	      $goView->graph->as_png("$EP::Config::h->{TMP}/$imageFile.png");
#	      $goView->graph->as_imap("$EP::Config::h->{TMP}/$imageFile.map");
	      $imageFile = $goView->showGraph;
	      $HTML_OUTPUT .= $goView->imageMap;
	  }
      }
      
      $HTML_OUTPUT .= table ( { -border => 0, -align => "left" },
			      Tr ( { -align => "center" },
				   td ( "legend: " ),
				   td ( { -bgcolor => "red" },
					  "<=1e-10"),
				     td( { -bgcolor => "orange"},
					 "1e-10 to 1e-8"),
				     td( { -bgcolor => "yellow"},
					 "1e-8 to 1e-6"), 
				     td( { -bgcolor => "#00ff33"},
					 "1e-6 to 1e-4"), 
				     td( { -bgcolor => "#7fffff"},
					 "1e-4 to 0.01"), 
				     td( { -bgcolor => "#5099ff"},
					 "0.01 to 0.05"), 
				     td( { -bgcolor => "#ebc79e"},
					 ">0.05")
				   )
			      );
      $HTML_OUTPUT .= "\n<br><br>\n";

      $HTML_OUTPUT .= h3( { -align => "center" },
			  "Most significant GO categories",
			  );
      $HTML_OUTPUT .= "<table border=1 cellpadding=0><tr><th>Gene Ontology term</th><th>Cluster frequence</th><th>Genome frequence</th><th>P-value</th><th>Genes annotated</th></tr>";
      # print list with go terms and associated genes
      foreach my $pvalue (@pvalues) {
	  next if $pvalue->{CORRECTED_PVALUE} > $cutoff;

	  # create html formatted links for nodes and genes
	  my $goId = $pvalue->{NODE}->goid;
	  my $nodeUrl = $EP::Config::h->{SPECIES_METADATA}->{$self->{species}}->{nodeUrl};
	  $nodeUrl =~ s/<REPLACE_THIS>/$goId/;
	  
	  foreach my $gene (keys %{$pvalue->{ANNOTATED_GENES}} ) {
	      my $StdGene = $self->{annotationProvider}->standardNameByName( $pvalue->{ANNOTATED_GENES}{$gene} );
	      my $geneUrl = $EP::Config::h->{SPECIES_METADATA}->{$self->{species}}->{geneUrl};
	      $geneUrl =~ s/<REPLACE_THIS>/$pvalue->{ANNOTATED_GENES}{$gene}/;
	      $pvalue->{ANNOTATED_GENES}->{$gene} = "<a href=$geneUrl TARGET=_NEW>$pvalue->{ANNOTATED_GENES}{$gene}</a>";
	  }
	  
	  $HTML_OUTPUT .= "<tr><td><a href=$nodeUrl TARGET=_NEW>" . $pvalue->{NODE}->term . "</a></td>";
	  $HTML_OUTPUT .= "<td>" . $pvalue->{NUM_ANNOTATIONS} . " out of " . scalar(@list) . " genes</td>";
	  $HTML_OUTPUT .= "<td>" . $pvalue->{TOTAL_NUM_ANNOTATIONS} . " out of " . $self->{termFinder}->__totalNumGenes  . " annotated genes</td>";
	  $HTML_OUTPUT .= "<td>" . sprintf( "%.3g", $pvalue->{CORRECTED_PVALUE} ) . "</td>";
	  $HTML_OUTPUT .= "<td>" . join( " ", values( %{$pvalue->{ANNOTATED_GENES}} ) ) . "</td></tr>";
      }
      $HTML_OUTPUT .= "</table>";
  

  # unused gene names
  if (@ambiguous) {
    $HTML_OUTPUT .= "<p>The following gene(s) are ambiguously named, and so will not be used:<br>\n";
    $HTML_OUTPUT .= join(" ", @ambiguous) . "</p>\n";
  }
  if (@notFound) {
    $HTML_OUTPUT .= "<p>The following gene(s) were not recognized, and will not be considered:<br>\n";
    $HTML_OUTPUT .= join(" ", @notFound) . "</p>\n";
  }

  $HTML_OUTPUT .= end_html();

  return $HTML_OUTPUT;
}

sub svg_html {
  my ($svg_file, $width, $height) = @_;
  my $svg_html = <<__SVG_APPLET__;
    <EMBED ALIGN="center" SRC="$svg_file" NAME="SVGEmbed" HEIGHT="$height" WIDTH="$width" TYPE="image/svg-xml" PLUGINSPAGE="http://www.adobe.com/svg/viewer/install/">
__SVG_APPLET__

  return $svg_html;
}

1;
